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ABSTRACT 

The observed gas-phase molecular inventory of hot cores is believed to be 
significantly impacted by the products of chemistry in interstellar ices. In this 
study, we report the construction of a full macroscopic Monte Carlo model of 
both the gas-phase chemistry and the chemistry occurring in the icy mantles 
of interstellar grains. Our model treats icy grain mantles in a layer-by-layer 
manner, which incorporates laboratory data on ice desorption correctly. The ice 
treatment includes a distinction between a reactive ice surface and an inert bulk. 
The treatment also distinguishes between zeroth and first order desorption, and 
includes the entrapment of volatile species in more refractory ice mantles. We 
apply the model to the investigation of the chemistry in hot cores, in which a 
thick ice mantle built up during the previous cold phase of protostellar evolution 
undergoes surface reactions and is eventually evaporated. For the first time, 
the impact of a detailed multilayer approach to grain mantle formation on the 
warm-up chemistry is explored. The use of a multilayer ice structure has a mixed 
impact on the abundances of organic species formed during the warm-up phase. 
For example, the abundance of gaseous HCOOCH 3 is lower in the multilayer 
model than in previous grain models that do not distinguish between layers (so- 
called "two phase" models). Other gaseous organic species formed in the warm- 
up phase are affected slightly. Finally, we find that the entrapment of volatile 
species in water ice can explain the two-jump behavior of H 2 CO previously found 
in observations of protostars. 
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Subject headings: astrochemistry, stars: formation, ISM: clouds, ISM: abun- 
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Introduction 



Interstellar molecules play an important role in the interstellar me dium (ISM) and exist 

both in the interstellar gas and in ice mantles on interstellar grains (e.g.. lHerbst &: van Dishoeckl 

2009). While many observed molecules are formed via gas-phase che mical reactions (IHerbst fc Klemperer 
19731 ). there are important exceptions such as molecular hydrogen (IGould fc Salpeterlll963l ) 
and complex organic m olecules widely observed in warm interstellar environments (e.g., 



Garrod &: Herbst 



20061 ). Even abundances of simple gas-phase species such as N 2 H + , HCO 



and C O in some cases may be controlled by chemistry in interstellar ices (jVasyunina et al. 
20121 ). 



Interstellar ices have been investigated extensively during the last four decades since 
the discovery of the 3 .1 /zm solid water infrared absorption band in Orion BN/KL by 



1984 1. CH 3 OH (lOrim et al 



ices to be CO ( 


Lacv et al. 


(de Graauw et al. 


1996). 



NH 3 , and H 2 CO (Ivan Dishoeckl |2004| ). The most abundant ice compound in cold sources is 
H 2 0, with a fractional abundance with respect to H 2 of ~10~ 4 , which is comparable to the 
concentration of CO, the most abundant gas-phase molecule. The next most abundant ice 
compounds are CO and C0 2 , with abundances relative to solid water of 0.1 — 0.3. Other 
species such as NH3, CH 4, and H 2 CO have abundances up to a few percent with respect to 
water (IGibb et al.ll2004j ). A comparison of ice inventories around high-mass and low-mass 
protostars shows a qualitative similarity of the ice composition in both types of objects. 
Quantitatively, however, the ice composition is somewhat different. In particular, the frac- 
tion of CO and C O? in the ices aro und high-mass protostars is lower than in ices around 
low- mass objects (jOberg et al.ll201ll ). possibly reflecting the volatility of these two species. 



Analysis of the shapes of absorption bands also yields constraints on the structure of 
interstellar ices and their evolution during the development of a protostar. Thick ice mantles 
typi cally consist of two major components: a water-rich polar ice and a water-poor apolar 
ice (ITielens et al.llf99ll ). While C0 2 normally resides in the polar component, the CO ice lies 
in the polar water-poor component. At late stages of star formation, the C0 2 :H 2 mixture 
is segregated forming C0 2 "island s" within the water ice, while the apolar CO component 
is gradually desorbed into the gas (IPontoppidan et al.ll2008l ; lOberg et al.ll201ll ). 
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Ices manifest themselves in observations of gas-phase molecules in vari ous star-forming 



enviro nme nts. Observations of Class and Class I protostars performed by iJprgensen et al. 



( 120021 ) and lSchoier et al.l (120021 ) indicate that the gas-phase abundance of the simple diatomic 
molecule CO might require entr apment of CO mole cules in w ater ice or othe r relat ively 
refractory compon ents of the ice (ICuppen et al.ll201ll ). Work by ICeccarelli et al.l (120011 ) and 
Maret et al.l (120041 ) show that formaldehyde H2CO might be trapped in the ice as well. In 
hot cores, which represent a l ater stage of protoste l lar evolution, complex org anic molecules 
are commonly observed (e.g.. iBottinelli et al.ll2004l . 120071 ; iBisschop et al.ll20071 ). The organic 
chemistry in hot cores is driven by heating and com plete evaporation o f ices formed during 



the previous cold stage of prestellar evolution (e.g., iBrown et al.l Il988l ; ICaselli et al.l 11993 
Garrod & Herbstl 120061 ) . Some organic molecules, such as methanol, are formed primarily 
on ices and then evaporate into the gas. For other molecules such as formic acid or methyl 
formate, surface formation routes are assist ed by gas-phase chemis try, which is activated by 
chemical precursors ejected from ices (e.g., iGarrod fc Herbstl 120061 ). 



The observational data on ices has given impetus to laboratory studies on interstellar 
ice analogs. Many studies done during the last decade were directed towards an under- 



standing of t he dynamics of ice mixtures in 



experiments (jCollings et al.l l2004j; lOberg et a 



PD ( temp erature programmed desorption) 



20071 . 120091 ) and the interaction of ices with 



UV photons. Experiments by ICollings et al.l (120041 ) revealed different patterns of desorption 
for light volatile species (N 2 , O2, CO, CH4 etc.) and heavy species with high desorption 
energies (H 2 0, NH 3 , CH3OH). In general, the desorption behavior of species in a mixed ice 
is complex, and depends on ice composition and thickness, as well as on the heating rate and 
any phase change in the ice during heating. With a high degree of simplification, one can say 
that light species can be trapped within more refractory ice, and exhibit at least a two-stage 
desorption consisting of desorption from the ice surface at their sublimation temperature, and 
co-desorption of entrapped volatile species with more refractory ice. Refractory species ex- 
hibit a much simpler desorption pattern, corresponding to their temperature of evaporation. 
In general, the desorption behavior of species from a thick ice typically shows zeroth-order 
behavior, in which the rate of desorption doesn't depe nd on total amount of ice molecules, 
while for thin ices, desorption has a first-order nature (jSomorjai fc Bentlll995l ). 



The majority of t heoretical astrochemic al models utili ze the treatment of grain s urface 
chemistry proposed in lHasegawa et al.l (119921 ) and earlier in iPickles fc Williams! ( 119771 ) . This 
is a very simplified approach in which no details of ice structure and physics is included. 
The numerical description of grain surface chemistry in such models is based on chemi- 
cal rate equations, which were later proven to have limited applicability to grain-surface 
chemistry, a problem that can be solved by using computer-intensive stochastic approaches 
( IVasyunin et al.l 120091 ). There are only a limited number of astrochemical studies where 
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more detailed models of grain surface chemistry have been used. In lHasegawa fe Herbst 
( 119931 ). rate equations based on a three-phase model consisting of g as, granular surfa c e, and 
inert bulk were introduced, but rarely used in subsequent studies. iGarrod fe Paulyl (120111 ) 
presented an updated version of this model, which includes the treatment of stocha stic pro- 
cesses (see below) via a semi-empirical mod ifi ed rate-equat i on ap proach proposed by IGarrod 
( 120081 ) and validated in IGarrod et al.l (120091 ) . iTaquet et al.l ( 120121 ) presented a similar model, 
in which the porosity of interstellar ices is also taken into account, and traditional rate equa- 
tions are utilized. The important feature of the three-phase models proposed so far is that, 
to the best of our knowledge, they are capable of building up an ice mantle in a layer-by-layer 
manner, but unable to evaporate it in the same manner, because only abundances of the 
outermost monolayer are considered. Below this layer, such models only remember average 
bulk abundances, without distinction between monolayers. The only model of diffusion and 



desorption of species trapped in a refractory ice is presented in a study of iFayolle et al. 
( 120111 ). but it does not include any chemistry. 



The first serious stochastic models were formulated by ICharnleyl (119981 . 120011 ). who 
introduced Monte Car l o mod els of grain surface chemistry based on a stochastic simulation 
algorithm by I Gillespie! (Il976l ). In these two papers, gas phase and grain surface chemistries 
were considered separately. The first Monte Carlo model of coupled gas-grain chemistry 



based on a large astrochemical network was presented by IVasyunin et al.l (120091 ). In this 



model, grain surface chemistry was also considered without a distinction between the reactive 
surface and an inert bulk made of multiple monolayers. 

A complex treatment of both the structure and molecular environment in all layers 
of grain mantles is embodied in models b ased on a microscopic Monte Carlo approach 
(jChang et al.l 120051 ; ICuppen fe Herbstl 120071 ). Here, the motions of individual atoms and 
molecules on grain ices are modeled. This type of model can include detailed physics of a 
grain mantle: layers, porosity, bulk diffusion etc. The main limitation of these models is 
computational, which renders the incorporation of large chemical networks and simulations 
of chemistry over long periods of time (more than 10 5 year) difficult. 

The aim of this study is to investigate the chemistry in hot cores affected by evaporation 
of ices formed during the previous cold stages of star formation with a new model, which 
includes recent d iscoveries in the ice ph ysics. For this purpose, a macroscopic Monte Carlo 
code reported in IVasyunin et al.l ( 120091 ) has been modified. In addition to the correct treat- 
ment of stochastic effects and a full gas-grain chemistry, it now includes several important 
properties of ices: layered structure, entrapment of volatile species in refractory water ice, 
and layer-by-layer desorption. With this model, we follow the chemistry starting from ice 
build-up during the cold collapse phase from semi-diffuse to dense cloud followed by the 



-5 - 



warm-up stage as a hot core develops. 

The outline of the paper is as follows. In Section &\ we describe our chemical and 
physical models, the algorithm of the model, and its numerical tests. In section [31 modeling 
results are presented. Section H]is devoted to a comparison of results with observational data 
and a discussion. A summary is given in section |5j Finally, the mathematical details of the 
algorithm used are presented in an Appendix. 



in 
In 



2. Modeling 

2.1. Modification of the Monte Carlo algorithm: the MONACO code 

In this study, we uti lize an extended version of the Monte Carlo approac h first describe d 



Vasvunin et al. 



Vasyunin et al. 



(2009), which is based on the classical Gillespie algorithm (lGillespidll976l ). 
(120091 ). the Monte Carlo algorithm was applied to a two- phase astrochem- 
ical system consisting of the gas phase and an undifferentiated grain mantle. The fact that a 
real grain mantle may be thick and consist of multiple monolayers, with just the outermost 
monolayer (or a few outer layers) chemically reactive while the bulk of grain mantle remains 
inert, was not taken into account. Here, we extend our algorithm in order to include a more 
realistic representation of a grain mantle in the model. 

Since our approach is macroscopic, we cannot incorporate many realistic microscopic 
details of the structure of an ice mantle such as surface roughness or bulk porosity, which 
can be taken into account in microscopic Monte Carlo simulations. Instead, we consider 
a very simplified model, although one that is more detailed than previous macroscopic ap- 
proaches. In our macroscopic Monte Carlo model, molecules are just logical entities that 
obey certain rules. We assume that: a) an ice mantle can consist of multiple monolayers; 
b) the formation of the (M+l)th monolayer starts only after the Mth monolayer is fully 
formed; c) only a certain number of outer monolayers are chemically active; i.e. chemical 
reactions and desorption are only possible between molecules belonging to outer monolayers. 
The exact number of chemically active monolayers is a subject of calibration of the model 
vs. laboratory data, as discussed in Section 12. 2| d) the number of molecules on the grain 
surface can be increased by accretion from the gas, decreased by desorption into the gas, and 
either increased or decreased by a surface chemical reaction depending upon the numbers of 
reactants or products; e) there is no porosity or internal motions of molecules in the bulk 
below the chemically active monolayers; f ) only molecules belonging to the chemically active 
monolayers can desorb. A detailed description of the algorithm is given in the Appendix. 



The important advantage of the developed algorithm is that it not only allows us to 



-6 - 



grow the mantle in a layer-by- layer manner, but makes it possible for it to be evaporated 
in a layer-by-layer manner as well, contrary to previous three-phase models. This feature 
is important for studies of chemistry during the warm-up phase; e.g., in hot corinos and 
cores. Because of this feature, we call our model "multiphase" in order to make a distinction 
between our model and previous three-phase models. 

The described algorithm was implemented in the MONACO code (acronym for "the 
MONte cArlo COde"), which is written in Fortran 90. The code is capable of simulating 
gas-grain chemistry under time-dependent physical conditions. Although the code is com- 
putationally expensive, it is efficient enough to be run on desktop workstations, albeit for 
large periods of time. 



2.2. Behavior of the model vs. laboratory data on ices 

It is important to benchmark the proposed algorithm against the results of laboratory 
studies of ices. The treatment of desorption and entrapment of species in the MONACO 
algorithm is of special interest, since this part of the algorithm is comp letely new. For 



this purpose, we used experimental results obtained by iFayolle et al.l ( 120111 ). In that study, 
TPD experiments involving well-mixed binary and tertiary ice mixtures (CO:H 2 0, C02:H 2 0, 
CO:C02:H 2 0) were performed with the CRYOPAD experimental setup. Experiments mea- 
sured the entrapment efficiency for the volatiles CO and CO2 within water ice. The entrap- 
ment efficiency is here the fraction of CO or CO2 with respect to its initial content in ice 
mixture, which remains trapped within the water ice, after the ice mixture is heated above 
the CO or CO2 sublimation temperature, but below the water sublimation temperature. 

It was found that CO and CO2 each have two desorption peaks, one of which is caused by 
desorption of these volatile species from surface layers of ice, and the other by co-desorption 
of refractory water ice and entrapped CO and CO2 molecules. The first peaks occur via 
first-order desorption, while the second peaks indicate zeroth-order kinetics, in which des- 
orbed molecules are replaced by species from within the bulk. The efficiency of entrapment 
is higher for CO2 than for CO. In general, entrapment is more efficient in thick ices (several 
tens of monolayers) than in thin ice films (6-16 monolayers) . A higher fraction of water 



molecules in ice also increases the efficiency of entrapment. IFayolle et al.l ( 120111 ) also de- 
duced that CO2 can diffuse out of the ice bulk to the surface and then evaporate only from 
several outermost monolayers of bulk ice, and there is a sharp boundary between monolayers 
available for diffusion and the rest of the bulk. Such a complex behavior cannot be fully 
reproduced within the framework of the simple MONACO algorithm, which totally neglects 
diffusion in the bulk ice. To understand the differences between our model and the lab 
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data, we simulated TPD experiment No. 19 from iFayolle et al.l (120111 ). in which desorption 
vs temperature was studied in a 30 ML (monolayer) tertiary ice mixture of H 2 0:C02:CO 
(20:1:1). In addition, we compared the results obtained wit h the MONACO code with an 
extended three-phase model proposed in lFayolle et al.l ( 1201 ll ) for the case of a thick 100 ML 
binary ice with a C02:H20 ratio of 1:5. There are no experimental data for such a thick ice 
in that pape r, but the ext e nded three-phase model is based on the whole set of experiments 
described in lFayolle et al.l ( 120111 ). and we believe that it gives a reasonable est i matio n of en- 
trapment efficiency for thick ice. For the case of experiment 19, IFayolle et al.l ( 120111 ) deduce 
an entrapment efficiency of 84%, for either the CO or CO2 species, while for the case of thick 
ice their model gives an entrapment efficiency of 95% for CO2. The second case should also 
be considered as astrochemically relevant despite its lack of CO, because the thickness of 
the interstellar ices is estimated to be ~100 ML and the relative fractio n of CO and CO2 to 
water varies between 1:10 and 1:3 ( IGibb et al.l 12004 : lOberg et al.ll201ll ). 



To model the TPD experiments, we set the initial ice com position in the MO NACO 



model according to those used in corresponding experiments of IFayolle et al.l (120111 ). The 
temperature in the model is then linearly increased over the timescale used in the exper- 
iments. In contrast to actual simulations described later in this study, no chemical reac- 
tions are considered. Therefore, the desorption e nergies of the species ar e the only param- 
eters needed. These ene rgies are 1150 K for CO ( G arrod &: Herbstl 120061 ) . 2850 K for CO2 



dGarrod fc Herbstl 120061 ) and 5700 K for H 2 dFraser et al.ll200lh . Figure □ represents the 
results of simulations of the TPD experiment with the tertiary C02:CO:H 2 ice mixture 
performed with the MONACO code in which one outermost monolayer is chemically active 
(left panel), or four outer monolayers are chemically activ e (right panel). From the com- 
parison of the figure with Figure 4a in IFayolle et al.l (1201 ll ). one can see that the results of 
the numerical TPD experiment obtained with the MONACO code are qualitatively similar 
to the real experimental results. Both CO (solid black) and C0 2 (dashed red) exhibit two 
peaks during the ice warm-up. The first peaks correspond to the evaporation of these volatile 
species from the ice surface, while the second peaks are caused by the co-desorption of en- 
trapped volatiles and refractory water ice. The shapes of the surface desorption peaks and 
bulk desorption peaks confirm that the MONACO code treats zeroth- and first-order desorp- 
tion correctly. The desorpti on peak of water is located in our model at higher temperatures 
than in the exp eriments of IFayolle et al.l (120111 ) because the desorption energ y of water we 



adopt following iFraser et al.l (120011 ) is somewhat higher than those found by IFayolle et al. 
(120 111 ) in their experiments. 



Let us now consider the quantitative agreement of laboratory TPD data and TPD 
simulations with our code. In Figure [2], the simulated fraction of CO2 in the 100 ML 
C02:H20 ice mixture versus temperature is shown vs temperature. Once again, the left 
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panel corresponds to our simulation with one chemically active monolayer, while the right 
panel represents results obtained with four chemically active monolayers. Using this plot, it 
is possible to estimate the entrapment efficiency of CO2 in the simulated ice sample. For that, 
one should divide the value of C02:H 2 ratio in the temperature range 60-100 K by the value 
of C02:H20 fracti on at tempera t ures b elow 50 K. A comparison of the left panel of Figure [2] 
with Figure 8 in iFayolle et al.l (120111 ) reveals that in our simulation with one outermost 
chemically active monolayer, the entrapme nt efficiency of volati le CO2 is significantly higher 
than in the extended three-phase model of IFayolle et al.l (1201 If ) ■ This fact can also be seen 
when comparing the left upper panel of Figure[Tjwith Figure 4a in lFayolle et al.l ( 1201 ll ). where 
the relative heights of the surface desorption peaks and co-desorption peaks are smaller in 
the simulation than in the real TPD experiment. Specifically, the entrapment efficiency in 
the simulation of the 30 ML tertiary ice COiCC^H^O is 96.3% vs 84% in the experiment 
and 98.8% vs 95% for the 100 ML binary C02:H 2 ice. We conclude that in the absence of 
bulk diffusion only a very small fraction of volatile species is available for desorption from 
the ice when the temperature is below the water evaporation threshold. 

To overcome this problem and make our modeling results closer to the laboratory data, 
we changed our initial assumption that only the outermost monolayer of the ice mantle 
is chemically active. It was found that by making the four upper monolayers of ice man- 
tle chemically active, we can make our modeling results much closer to the results from 
Fayolle et al.l ( 120111 ) for the astrophysically important case of the 100 ML binary C02:H 2 
ice. This assumption also helps to reproduce the TPD experiment with the 30 ML ter- 
tiary ice. The results of our simulations with the tuned-up MONACO code are presented 
in the right panels of Figs. [TJ and [2], and may be compared with Figure 4a and Figure 8 



111 



Fayolle et al.l ( 120111 ). Our revised values for the entrapment efficiency are 85% for the 



tertiary ice and 95.2% for the binary ice, which lie very close to experimental and model 
values. It is important to note that our revision of the code cannot fully compensate for the 
absence of bulk diffusion. In general, one can expect that the revised code will overestimate 
the efficiency for the entrapment of volat ile species in t he ca se of ices with a low fraction 
of refractory water (see, e.g., Table 1 in IFayolle et al.l ( 1201 ll )). However, for the purpose 
of this study, in which the evaporation of thick ice with water as the major component is 
investigated, the adopted revision is reasonable. 
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2.3. Chemical model 



In this study, we utilize a full gas-grain chemical m odel, in which the gas- grain chemical 
network is a version of the KIDA database published in ISemenov et al.l ( l2010f rl with several 
minor updates concerning missing desorption processes for several molecules in the original 
file. The KIDA network in turn, is a successor of the OSU astrochemical database, which can 



be found at |http://ww w. physics.ohio-state.edu/~eric/research.html (IWakelam et al.ll2012l ). 
In total, the chemical network includes 662 species and 5693 chemical processes. The surface 
network does no t include the many additional complex molecules in the most recent OSU 
network used by lGarrod et al.l ( 120081 ) . Nevertheless, it is sufficient for this work, since in the 
warm-up phase we limit our study of complex organic species to dimethyl ether (CH3OCH3), 
formic acid (HCOOH), and methyl formate (HCOOCH3). 

There are two major modifications introduced to the original rate file. First, the disso- 
ciative recombin ation channels o f proto nated dimethyl ether were replaced with those taken 
from Table 4 of lHamberg et al.l ( |2010| ). Second, the formation of protonated methyl for- 
m ate via the gas-ph ase route CH 5 + + H 2 CO — > H 5 C 2 0^ + H 2 was disabled according 
to iHorn et al.l ( 120041 ) . Also, the desorpt ion energy of molecular hydrogen was chosen to 



be 314 K according to iKatz et al.l ( 119991 ) . The gas-phase and grain surface chemistry are 
connected by accretion and desorption processes. For accretion, the sticking probability for 
all neutral species except for H, H 2 and He is unity. Molecular hydrogen, helium and ions 
are not allowed to stick on the grain surface, although selected ion reactions with negatively 
charged grains are included. The non-sticking of molecular hydrogen in our model is dictated 
by computational problems. Including the accretion of molecular hydrogen would slow down 
our Monte Carlo simulations by a factor of ~10000, and make them unfeasible. The role of 



reacti ons involving radicals and H 2 on surfaces has been investigated by lHasegawa fe Herbst 



( 119931 ) in a model that uses simple rate equations with a common barrier for all such reac- 
tions. In general, hydrogenation on surfaces is handled efficiently by association reactions 
with atomic hydrogen. In our network, we produce water via the well-know n sequence of 
reactions between atomic oxygen and atomic hydrogen (e.g.. lDulieu et al.ll2010l ). As is shown 
later in this study, the water on grain surfaces is produced in abundance by this mechanism 
at the studied temperatures. Desorption processes enabled in our model include thermal 
evaporation, cosmic ray induced desorption, and photodesorption. Given the uncertainty 
of laboratory experiments and lack of experimental data for the majority of molecules, we 
adopt ed a photodesorption yield of 10~ 3 for all surface species according to lOberg et al. 



( 120071 ) . Only molecules in chemically active monolayers are subject to desorption. 



1 kida.obs.u-bordeauxl.fr /uploads /models /benchmark_20 1 .dat 
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We adopted a cosmic ray ionization rate for H2 of £=1.3-10" 



17 s- 1 . 



The single granular 



radius in our mode l is 10 5 cm. The elemental abundances correspond to set EAl from 



Wakelam fe Herbstl (120081 ). All elements initially are in atomic form except for hydrogen, 
which is initially bound up in H 2 . 

We consider only diffusive surface chemistry via the Langmuir-Hins helwood mecha 



nism. The surface mobility of species is described using the formalism from lHasegawa et al. 



(119921 ). According to this formalism, there are two sources of mobility for species on grain 
surfaces: thermal hopping and quantum tunneling through diffusion barriers for the lightest 
species — H and H2. The main control parameters for the diffusive chemistry — the diffu- 
sion/desorption energy ratio E^/Ed and the thickness of the potential barrie r for diffusion - 



are k nown poo rly. In the literature, ad opted values for E h /E D vary from 0.3 ( lHasegawa et al. 



19921 ) to 77 (IRuffle fe Her 



DSt 



20001 ). Tunneling for the lightest species is assumed to b e 



efficient in Hasegawa et al. ( 19921 ). but late r studies do not confirm that (Katz et al. 1999 ) 



although some recent models also include it (ICazaux fe Tielensll2004l ; llqbal et al.ll2012l ). Hav- 
ing in mind these uncertainties, we utilized three sets of control parameters for the surface 
mobility of species in this study. In the first set, we assume a diffusion/desorption energy 
ratio of 0.3 and efficient tunneling for H and H 2 on surfaces through thin rectangular poten- 
tial barriers of width lA. In the second and third sets, tunneling for light species is switched 
off, with Eb/E-o=0.5 in the second set, and 0.77 in the third s et. Therefore, the first and 
the third sets co r respon d to extreme values considered before in lHasegawa et al.l (119921 ) and 
Ruffle fe Herbstl (120001) . The second set re presents an intermediate case and was utilized 



previously in, e.g., iGarrod fe Herbstl (120061 ) . The three sets of control parameters for the 



mobility of surface species are identically implemented in the MONACO code and in a stan- 
dard two-phase rate-equation code (RE), which is used for comparison of the new multilayer 
approach with traditional models. 

For all three sets of control parameters, the rates of all surface reactions that have 
activation barriers were calculated assuming quantum tunneling throu g h rec tangular acti- 
vation barriers of width 1 A according to Eq. (6) in lHasegawa et al.l ( 1992). The set of 



activa tion barriers of surface reactions originates from the earlier work of IGarrod fe Herbst 



(2006), whi ch in turn is pa rt ially based on the compilation of sets from lHasegawa et al 
( 119921 ) and IRuffle fe Herbstl (|2000| ). The majority of activation barriers are taken from 



analogous gas-phase reactions. This fact motivated our decision not to include a so-called 
reaction-diffusion competition for surface reactions with activation barriers in the model 



( iHerbst fe Millar! 120081 ). Competition only pertains to diffusive surface reactions, not to 
gas-phase ones. Therefore, the values of activation barriers inherited from measurements 
of gas-phase reactions may not be valid when reaction-diffusion competition is enabled for 
surface reactions. Instead, barriers obtained directly for surface reactions should be used 
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Cuppen et al.ll2009l ; IChang &: Herbstll2012l ). These data are available only for several 



(e.g., 

reactions, and their incorporation into a full-scale network of reactions is problematic. 



We do not use the modified rate equations proposed by iGarrodl (120081 ) in this study 
despite their relative success in other studies, because they were formulated only for those 
two-body grain surface reactions, for which reactants appear on a grain surface independently 
of each other. If a surface reaction has products that are reactants for other surface reactions, 
such as photodissociation, the crucial condition of independence is violated. Further study 
on this point is warranted. 



2.3.1. Multiple layers of ice mantles 



According to lHasegawa et al.l (119921 ). the mobility of species on grain surfaces can be 
converted to the rate of a surface two-body reaction between species i and j using the 
expression: 

ru = kMjNs + rijN^im)^) (l) 

where rh op is the hopping rate or tunneling rate for a species over or through the diffusion 
barrier to a nearest neighbor site, k%j is a factor < 1 related to tunneling under any activation 
energy of the reaction, if it exists, (rii) and (rij) are the average abundances of species i and 
j on a grain surface, and N s is the total number of adsorption sites on the surface. 

This expression is reasonable when the total number of species on a grain surface N tot 
does not exceed iV s . However, if the total number of species on a grain exceeds N s , some 
fraction of molecules in the ice lies buried within the mantle and might not be available for 
diffusion. In the case of a thick mantle, the majority of species are buried. Hence, if we 
assume that only molecules of one surface monolayer are reactive, we must correct expression 
(ffl) to 

rf w = rij • [min(N s /N tot , l)] 2 . (2) 

For the case of unimolecular surface photoreactions, if we believe that they occur only in the 
outermost monolayer, expression ([2]) should be modified to 

r% ew = r ir min(N s /N tot ,l) (3) 

These corrections, in principle, should be implemented in the two-phase rate-equation model. 
However, we did not include it in in our RE model so that our two-phase model could resemble 
as closely as possible those used in previous studies. 



Since we assume four upper monolayers to be chemically active according to Section I2T2| 
but also wish to investigate the impact of corrections ([2]) and ([3]) on the modeling results, 
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we consider two variations of our MONACO model: MC-1 in which corrections (|2]) and 
([3]) are implemented, and MC-4, in which we omit this correction. In Model MC-1, the 
corrections lead to the result that when the mantle thickness is four monolayers or more, 
the rates of unimolecular and two-body surface reactions are slower by a factor of 4 and 
16, respectively, than in model MC-4. Thus, we have three models - RE, MC-1 and MC-4 
- with physically distinct assumptions, which allow us to investigate the differences among 
assorted two-phase and multiphase models, each of which was run with the three sets of 
E h /E D discussed previously. 



2.4. Physical model 



For this study, we adopt a schematic two-stage model of formation of a hot core aiming 
to mimic the evolution of a gas parcel from the diffuse cloud phase to the hot core phase. The 
idea of the m odel is similar to those described in lGarrod fc Herbstl ( 120061 ). and inherited from 
earlier work ( jViti et al.l 120041). The model consists of two stages. In the first "cold" stage, a 
free-fall collapse occurs (jSpitzerl Il978l ; iBrown et al.l 119881 ). Physically, we attribute the cold 
collapse phase in our model to the formation of a cold prestellar cloud from a translucent 
cloud. The collapse starts at a gas density nn = 3 x 10 3 cm -3 and a visual extinction of 
Ay=2 and continues until the desired final density is reached. The collapse continues until 
a gas density of 10 7 cm -3 is reached after 10 6 yr of evolution. Figure E] shows the change in 
the visual extinction, gas density, and temperature during the cold collapse. Despite the fact 
that free-fall collapse is isothermal, in our study we assume that the temperature linearly 
drops from 20 K to 10 K with time during the collapse. Physically, this corresponds to the 
fact that the efficiency of radiative heating of grains drops as the visual extinction increases. 
Given the number of uncertainties in our knowledge of the collapse and of the properties of 
interstellar grains, we believe that such a schematic approach is sufficient fo r our model, and 
a more detailed treatment of the dust temperature (IGarrod fc Pauryn201l[ ) is not necessary. 
The gas temperature is assumed always to be equal to the dust temperature in our model. 

The second stage of our model is the warm-up phase. Here we assume that as the gas 
and dust collapse to wards the protos tar, a quadratic warm-up fro m 10 K to 200 K over 
2 x 10 5 yr following IViti et al.l ( 120041 ) and IGarrod fc Herbstl ( 120061 ). which corresponds to 
intermediate mass star formation. In the warm-up phase, the density remains constant. 
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3. Results 

3.1. Cold Collapse phase: build-up of the ice mantle 

3.1.1. Early collapse stage 

The fractional abundances of selected simple gas-phase species during the whole cold 
collapse phase are shown in Figure H] for Model MC-1. This figure has three panels, one 
for each E^/Ed ratio considered. Note that there is little difference among the abundances 
in the three panels. The abundances are affected even less if we switch to model RE or 
MC-4, so nine panels are not required. Figure [5] contains fractional abundances for the 
major ice species as a function of time during the whole cold collapse stage. Here there 
are nine panels since differences among the two-phase rate-equation approach (RE) and the 
multilayer approaches (MC-1 and MC-4) can be seen. 

During the early stage of collapse (up to 2 x 10 5 yr), the chemical evolution of the major 
ice mantle species is qualitatively similar in the RE, MC-1 and MC-4 models (see each column 
in Figure |5]). This result is natural, because at this very early stage the thickness of the ice 
mantle is < 10 monolayers, as shown in Figure [61 and the number of molecules participating 
in surface chemistry in similar all models. On the other hand, the different E^/Ed ratios can 
cause more significant differences in the ice composition, although models with E^j Ed=0.3 
and 0.5 are quite similar during this early stage. The two major constituents of the ice 
in these models are carbon dioxide (CO2) and water (H 2 0). Indeed, more than 90% of 
the mantle consists of these two species at early times. Both molecules form on the grain 
surface in competition with each other via the reactions grOH + grCO — > grC0 2 + grH and 
grOH + grH — > grH 2 0. Here gr indicates that species are on the grain surface. Atomic 
hydrogen comes from the gas phase where it is produced via dissociation of H 2 by photons 
and cosmic rays (see Fig. [4]). About 50% of the carbon monoxide is also accreted from the 
gas, but it is also produced on grain surfaces via the routes grC + gr0 2 —> grCO + grO and 
grC + grOH — > grCO + grH. The common precursor of both C0 2 and H 2 0, hydroxyl radical 
OH, is mainly formed from the atomic species accreting from the gas. 

In the multiphase models with E^/E-q = 0.3, the competition between C0 2 and H 2 
formation is won by water from the earliest times until the end of cold collapse phase. 
This is the most abundant compound in the ice since the first built monolayer in both 
Monte Carlo models MC-1 and MC-4, as can be seen in Figure [5J left column and also 
Figures [7J and top row, middle panel. In the middle panels of these latter two figures, 
the molecular abundances are expressed in terms of fraction of the grain mantle, whereas in 
the left panels they are expressed in terms of fractions of individual monolayers. In the RE 



-14- 



model, CO2 is more abundant than H2O. We believe that this is due to stochastic effects in 
grain surface chemistr y, which cannot be treated correctly in models based on rate equations 



flVasvunin et all 120091 ) . 



In the models with E^/Ed = 0.5, the competition is won by CO2 only in the early 
stage. In these models, the quantum tunneling as a source of atomic hydrogen mobility 
on the surface is switched off, and the majority of H atoms desorb back to the gas before 
they react with OH. The heavier molecules OH and CO are better bound to the surface, 
and therefore they can react at T dust =17-20 K, the relevant early stage temperature range. 
Although an activation barrier of 80 K slows down this reaction, a combination of T~20 K, 
the relatively high mobility of CO at this temperature, and E^/Ed=0.5 makes this reaction 
the major route of CO2 formation. For the model with E^/E-q = 0.5 at T=20 K, the mobility 
of OH is negligible: ~10~ 19 sites s _1 , while the mobility of CO is ~1 site s _1 . 

The ice composition in models with E^/Ed = 0.77 (corresponding to low surface mo- 
bility) differs from other models even at the early stage, as can be seen in Figure [5j right 
column, as well as the bottom row of Figures [7] and El The largest constituent is carbon 
monoxide CO instead of water, while the fraction of CO2 in the entire mantle is always 
less than 10%. This dramatic difference is caused by the low mobility of surface species 
at E^/Ed = 0.77, which makes the grCO + grOH — > grC02 + grH reaction inefficient in 
contrast to the model with E^/ Ed=0.5. Thus, most of the CO remains in its pristine form 
in the ice mantle, where at the early stage it is formed in the grC + gr02 — > grCO + grO 
and grC + grOH — > grCO + grH reactions. 

Two other major species in the ice mantle, ammonia NH3 and methane CH4, have much 
smaller abundances than CO2 and H2O. Both NH3 and CH4 are formed via similar chains 
of reactions: grNH — > grNH 2 — > grNH 3 and grCH — > grCH 2 — > grCH 3 — > grCH 4 . Another 
abundant surface species at the early stage, molecular nitrogen (N 2 ) is mainly formed by the 
surface reactions grN + grN — > grN2 and grNH + grN — > grN2 + grH. Also, about 40% of N2 
accretes from the gas phase. 

Since visual extinction is low at early times (Ay < 3), the photodissociation of molecules 
on grain surfaces is quite active. At the same time, it does not strongly affect the chemical 
composition of the ice, because the major molecules are involved in formation-destruction 
cycles such as grH + grOH — > grH 2 0, grH 2 + hv — > grH + grOH. The net effect of such 
loops is mainly a somewhat delayed build-up of the mantle, which can be seen as a less steep 
slope in Figure El at early time. 
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3.1.2. Intermediate collapse stage 

After ~ 2 x 10 5 yr of collapse, differences between models with various E^/Ed ratios 
as well as between the two-phase rate equation-based and multiphase Monte Carlo-based 
models become more pronounced. The most noticeable transition is a sharp decrease in the 
fraction of carbon dioxide CO2. In models MC-1 and MC-4 with E^/E^ = 0.3, the 60 lowest 
monolayers of the mantle are rich in CO2, and the sharp fall-off of its fraction happens at 
t~ 5 x 10 5 yr, as can be seen in Figures [7] and |HJ For Monte Carlo models with E^/Ed = 0.5, 
only ~ 25 monolayers are rich in CO2, and the sharp fall-off of CO2 fraction happens earlier, 
at t~ 3 x 10 5 yr. The fall-off occurs as the reaction grOH + grCO — > grC02 + grH becomes 
inefficient, which occurs at 17 K if E^/Ed = 0.5 and 15 K if E^/En = 0.3. As the temperature 
falls, CO2 formation loses the competition with H2O formation. In models with E^/ Ed=0.3 
it happens later than in models with Ey,/ Ed=0.5 because of the lower surface mobility of 
species at higher E^/Ej). From this point of view, models with Eb/En=0.77, where CO2 
formation is not efficient from the very beginning, can be considered as models where surface 
mobility of CO is too low to form CO2 even at T~20 K. At the late stage, in models with 
E^/ Ed=0.3, the fraction of CO2 increases again because the abundance of atomic hydrogen 
drops in the gas phase with increasing density, and the OH on the surface is not consumed 
in H 2 formation efficiently. 

Interestingly, as can be seen on top left panel of the Figure in the RE model with 
£'b/£^D=0.3, CO2 remains the most abundant fraction of the mantle until the end of the 
simulation. This arises from the fact that in the two-phase RE model the entire mantle is 
accessible for surface photodissociation processes, contrary to the multiphase MC-1 and MC- 
4 models, where only the upper monolayers can be photodissociated. The significant fraction 
of newly formed water is dissociated back to H and OH even at the late stage, because of 
cosmic ray-induced photodissociation. OH, in turn, can react with abundant CO and form 
CO2, even though the OH+CO channel is relatively inefficient. In multiphase models MC-1 
and MC-4, the majority of water is buried within the thick mantle and is not accessible for 
photodissociation. This makes the final abundance of water higher in MC models than in the 
RE model and damps CO2 formation at the late stage. Since CO2 can be photodissociated 
too, in RE models with E^/ Ed=0.5 and Eb/Ej)=0.77 one can see a late decrease in the CO2 
abundance. There, the mobility of CO is not enough to resist photodissociation, and the 
majority of CO2 is converted back to CO. In the model with E^/ Ed=0.5, CO is ultimately 
converted to methanol, while in the model with E^j Ed=0.77 it remains in its original form. 

Another important chemical difference between the intermediate and early collapse 
stages in models MC-1 and MC-4 with Eb/Ei)=0.3 and Eb/E-[)=0.5 is the activation of 
the CO hydrogenation sequence grCO — > grHCO — > grH 2 CO — > grH 3 CO — > grCH 3 OH. 
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This very important process is not efficient at early time because the temperature is too 
high for H atoms to be efficient at surface hydrogenation and because the high abundance 
of OH destroys formaldehyde in the reaction grOH + grH 2 CO — > grHCO + grH 2 0. Once 
the temperature decreases to <15 K, the hydrogenation sequence becomes more effective. In 
models with E\> / ' Ed=Q .3 , CH3OH becomes abundant at ~ 5 x 10 5 yr, while its fraction in the 
newly formed mantle monolayers reaches ~30%. In models with E^/E-[)=0.5, on the other 
hand, the maximum abundance of methanol in these mantle monolayers is only ~10%, a 
number reached later, after ~ 7 x 10 5 yr of collapse. Formaldehyde (H 2 CO), the precursor of 
methanol, is also formed in ice in significant amounts in the intermediate and late stages. In 
the RE model with E^j Ed=0.3, the methanol fraction never becomes significant. Again, this 
is a manifestation of the two-phase approach, where methanol is efficiently photodissociated 
in the entire mantle to formaldehyde. About 25% of the H 2 CO, in turn, is photodissociated 
to CO and H 2 . Finally, CO is converted to C0 2 . In models with E h /E D =0. 77 methanol is 
not produced in appreciable amounts. 



3.1.3. Late collapse stage 

The late stage of collapse (t > 7 x 10 5 yr) is characterized by a rapid increase in the 
gas density and visual extinction Ay, as can be seen in Figure [3J In the gas phase, the 
abundances of species first reach values typical for cold prestellar clouds (see Figure M for 
t ~ 7 x 10 5 yr). Then, after 8 — 9 x 10 5 yr, gas-phase species rapidly deplete onto grains 
due to the high density and low temperature. Among these species is the most abundant 
gas phase molecule, carbon monoxide. The composition of ice monolayers formed during 
this period is characterized by an increasing fraction of CO, as seen in Figures [7| and [5J 
Correspondingly, the fraction of other ice constituents decreases in the outer monolayers. 
In all Monte Carlo models, CO becomes the most abundant ice species in the outer 10-30 
monolayers, depending on E^jE^. 

The analysis above shows that one can delineate three chemically different phases of 
the ice mantle build-up during the contraction of a cold core in multiphase models with 
E h /E D =0.3 and E h /E D =0.5. At t < 2 x 10 5 yr, the C0 2 -rich fraction of the ice mantle 
is formed, and is intimately mixed with water, which is formed at the same time. Thus, 
we can attribute t his fraction of the i ce to the polar C0 2 :H 2 component of ice inferred 



from observations (ITielens et al.lll99ll ; IQberg et al.ll201ll ). Between 2 x 10 5 yr and 8 x 10 
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yr, an H 2 0-rich ice is formed. Other major ice constituents at this time are CO, ammonia 
NH3, formaldehyde H 2 CO, and methanol CH3OH. The third phase of ice is formed during 
the late stage of collapse between 8 x 10 5 yr and 10 6 yr. This period is characterized by 
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the massive freeze-out of CO from the gas phase onto the grain surface. Thus, the outer 
layers of the ice mantle are dominated by solid CO and may be associated w ith the inferred 
apolar component of interstellar ices fjTielens et al.lll99ll ; lOberg et al.ll201ll ). This picture, 
generally consistent with observations, is not clearly pronounced in the two-phase RE models 
considered in this study. Instead, there are no inert ice monolayers the composition of which 
reflects the history of changing physical conditions during the ice formation. Rather, the 
entire mantle is available for chemical reactions all the time, which leads to a mixed ice 
composition. In models with E^/ Ed=0.77 the first C02-rich phase of the ice mantle is never 
produced. This failure to reproduce observational inference argues against an E^/Ed ratio 
equal to 0.77. 



3.1.4- Comparison with observations 
It is interesting to compare the modeling results with obse rvational data of co ld ices. 



To do this, we took the median ice composition from Table 2 of lOberg et all (120111 ). which 
contains sources involving low- and high-mass protostars as well as background stars. These 
data are based on the critical compilation of observational data obtained with the Spitzer 
and ISO infrared space telescopes, and recommended by the authors as a reference for ice 
formation models. For comparison with our simulations, we used the median ice composition 
for low-mass protostars because our temperature profile during collapse corresponds to the 
early stage of cloud contraction. Thus, there is no protostellar heating phase, which is 
probably responsible for the differences in ice composition between low-mass and high-mass 



protostars in later stages of evolution (lOberg et al.ll201ll ). The protostellar heating phase 



can rather be attributed to the warm-up stage of the model, where temperature increases 
from 10 K to 200 K occur. One should bear in mind that the observed ice composition in 
each category is obtained from the analysis of infrared absorption spectra, which are usually 
averaged over regions with different physical conditions. In our model, we consider just 
one parcel of a gas-grain mixture with homogeneous (although time- dependent) physical 
conditions. Therefore, some discrepancies between observed and modeled ice compositions 
are natural. 

Our time-dependent values for the fractional abundances of major ice species per mono- 
layer and per mantle are shown in Figures [7] and [S] for the MC-1 and MC-4 simulations 
respectively. We compare our results for overall abundances per mantle with observed me- 
dian abundances as a function of time to determine the optimum agreement time t C omp, when 
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the sum of squares of the residual differences is minimal: 

N 

tcomp : F fit = J2( X ° bs - Xr del (t)) 2 = min. (4) 
i=i 

Here, t comp is the time of best agreement, X° bs (t) is the observational abundance of species 
i with respect to the abundance of solid water, while X™° del (t) is the modeled abundance 
of the same species at time t. The smaller the value of F^ t , the closer the modeled ice 
composition to the observed composition. Histograms of the best fit to the ice composition 
for MC models are shown in the right column of Figures [7] and while histograms of the 
best fit to the observed ice composition for RE models are shown in Figure [10j 

The optimum time and ice composition corresponding to the minimum values of F^ t for 
different models are summarized in Table [1] along with the observed median values for low- 
mass protostars, with the water ice abundance fixed at 100. F&t changes with time smoothly, 
so one should consider the best fit time as a rough estimation, not a precise value. One can 
see that reasonable ice compositions are produced by MC and RE models with E^/ Ed=0.5. 
Models with E h / E D =0.3 underproduce solid CO. At the same time, the calculated amount of 
CO2 is very high in the RE model (170% of that of water), while in the MC models CH3OH 
is significantly overproduced. Specifically, in the MC E^/E^—0.3 models, the optimum CO2 
abundance lies between 22% and 32%, while the methanol abundance is ~34%-35%. In 
addition, the methane abundance in these models is 2-3 times higher than the observed 
value. Models with E h /E D =0.77, on the contrary, overproduce the solid CO fraction and 
greatly underproduce CO2 and methanol. Specifically, the CO fraction varies from ~50% in 
the RE model to ~50%-70% in the MC-1 and MC-4 models. At the same time, the C0 2 and 
CH3OH fractions are calculated to lie close to zero. The calculated methane abundance in 
models with E^/ Ev=0.77 is in reasonable agreement with observation, at only 1.5-2.0 times 
higher than the observed value. 

Of the three models that yield very low values of F^ t , the multiphase Monte Carlo 
models with Eb/En=0.5 produce mantle compositions closer to what is observed. Of these 
two, Model MC-1 produces the better fit as indicated by the lower value of F fit =0.03 for 
MC-1 vs. Ffi t =0.07 for MC-4. For MC-1, the modeled values are closer to median detected 
values, and deviate from them by a factor of 1.5-2.0, except for NH3. Model MC-4, on the 
other hand, overproduces methanol and underproduces CO2. Although the RE model with 
E h / E D =0.5 possesses the minimum value of F fit =0.02, it should be remembered that F fit 
characterizes the difference between modeled and observed abundances of ice constituents 
with respect to solid water. In the RE model with E\>/ Ed=0.5 an optimum ice composition is 
produced at an optimum time of ~ 4 x 10 5 yr, when the conditions have not yet reached high 
density and Ay, and the abundances of ice constituents with respect to hydrogen number 
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density are at least an order of magnitude too low to explain the observed values. 

Let us now compare the modeled abundances of selected gas-phase species shown in 
Figure with observational values. These abundances reach the peak values at ~6xl0 5 yr 
when the gas density reaches ~10 4 cm" 3 . Later, due to high density, the gas-phase species 
are depleted on grains rapidly, and their abundances get very low. Our modeled peak 
abundances of gas-phase species are in reasonable agreement with observational values for 
the dark cloud core TMC-1 (CP), as seen in Table [2J 



3.2. Warm-up phase of a hot core 

We now discuss the chemistry during the warm-up stage leading to the formation of a 
hot core, using both two-phase and multiphase models. As discussed in the previous section, 
we ran the RE, MC-1 and MC-4 models, each with the three different values of E^/E^. 

The chemical evolution during the warm-up phase is presented in Figures [TTJ [121 an d 
[TBI In our model, the warm-up stage begins after the cold collapse stage is finished; i.e., 
after 10 6 yr of evolution. In Figures HU [121 an d H3] the warm-up time starts from zero, 
but the zero point corresponds to the end of the cold collapse phase. The warm-up from 
10 K to 200 K occurs from yr to 2xl0 5 yr. Figure [TT] shows the evolution of gas-phase 
species that are also main constituents of ice mantles during the early warm-up stage, having 
accumulated on grain surfaces during the cold collapse phase. During the warm-up, they are 
eventually ejected into the gas phase by thermal and non-thermal mechanisms. In Figure [T2~| 
the evolution of gas-phase oxygen- containing organic species formed via a combination of 
gas-phase and grain-surface chemistry during the warm-up phase is shown. Figure [TBIexhibits 
the evolution of their grain surface counterparts. The noise in low abundance seen in the 
figures is an intrinsic feature of Monte Carlo simulations, which deal with whole numbers of 
species. In our model, an absolute abundance of unity corresponds to an abundance relative 
to the total number of hydrogen nuclei of ~10 -12 . 



3.2.1. Abundant species 
The top row of Figure [TTJ represents the results of the RE models . Our two-phase 



RE m odel with E\>/ En=0.5 essentially reproduces the results obtained by lGarrod & Herbst 



( 120061 ). who used a similar approach. All species except atomic hydrogen are desorbed into 
the gas phase from the ice mantle at their sublimation temperatures. In the RE model, 
species evaporate independently of each other, and evaporation happens in one step, after 
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which abundances of these rather stable species stay nearly constant until the end of the 
simulation. The abundances of the gas-phase species, once sublimation has occurred, are very 
similar to the abundances in the ice before sublimation occurs. The exception is ammonia, 
which drops by 3 orders of magnitude in less than 2 x 10 5 yr after sublimation in model RE 
with E h /E D =0.3. 

The middle and bottom rows of Figure [TT1 show the results obtained with the multiphase 
Monte Carlo models MC-1 and MC-4. The main difference from the RE model is the two- 
step picture of the sublimation of volatile species discussed earlier in this paper. The first 
sublimation step occurs when sublimation temperature of the volatile species is reached, 
at which time only the species from four chemically active monolayers desorb. The second 
step occurs at the sublimation temperature of water ice, which releases all trapped volatile 
species as well. The two-step phenomenon is a new phenomenon in warm-up chemistry 
(if not in the laboratory), completely missing in two-phase models. Interestingly, after 
the first sublimation step, CO2 in the gas phase has a similar abundance in models with 
E^/ Ed=0.77 despite the fact that there was no CO2 on the ice surface at the end of cold 
collapse phase. This CO2 is produced initially on the ice surface during the early warm-up 
when the dust temperature reaches 20 K. This early CO2 may contribute to the formation 
of a segregated layer of pure CO o observed towards protostars, and usually c onsidered as 
evidence of moderate ice heating ( IBoogert et al.ll2008t iPontoppidan et al.l 120081 ). 



3.2.2. Complex organic species 



In Figure H2J the gas-phase abundances of organic molecules produced mainly on warm 
ing ice surfaces followed by sublimation are presented as functions of time and rising temper 
atures. Methanol and formaldehyde, as precursors of these species, are also shown. Figure 
shows the time dependence of the abundances of these molecules in their solid phase, which 
are formed by radical-radical recombination during warm-up. The formati on of these species 



i n the framework of t wo-ph ase models has been d iscussed in detail by iGarrod fc Herbst 
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multilayer approach to surface chemistry affects the previously proposed mechanisms. 



According to the analysis performed by lGarrod fc Herbstl (120061 ). methyl formate HCOOCH3 
is formed primarily on grain surfaces via the reaction 



grHCO + grCH 3 -> grHCOOCH 3 , 
assisted by a gas-phase reaction: 

H 2 COH+ + H 2 CO -»■ H 2 COHOCH+ + hi/. 



(5) 



(6) 



In this work , we make no dist inction between the cis and trans forms of methyl formate, as 
discussed by lLaas et al.l ( )201ll ). Although formic acid HCOOH and dimethyl ether CH3OCH3 
also have surface formation routes: 



grHCO + grOH -> grHCOOH, (7) 

and 

grCH 3 + grCH 3 grCH 3 OCH 3 , (8) 

the majority of these two species are produced by gas-phase routes. During the early warm- 
up, HCOOH is formed in the gas phase via the dissociative recombination of HCOOH^ . The 
sharp rise in the HCOOH gas phase abundance during the late warm-up is caused by the 
gas-phase reaction 

OH + H 2 CO ->■ HCOOH + H, (9) 

which becomes effective after the hydroxyl radical desorbs from grains. For dimethyl ether, 
the main formation route is a gas-phase reaction of methanol and protonated methanol: 

CH3OH + CH 3 OH^ -»■ CH 3 OCH| + H 2 0, (10) 

followed by dissociative recombination. This reaction becomes important when methanol 
evaporates from the grain mantle. 

All models with Et,/En=0.77 are characterized by a low amount of formaldehyde and 
methanol in cold grain mantles as seen in panels of the right column of Figure [T3l As a result, 
the peak abundances of these species in the gas phase after sublimation during the warm-up 
phase are 1-4 orders of magnitude lower than in models with E^/Ed=0.3 and E^j Ed=0.5. 
Since the gas-phase route of formation for dimethyl ether (see eq. [TOl) becomes inefficient, the 
peak gas-phase fractional abundance of CH 3 OCH 3 in the warm-up phase drops by 5 orders 
of magnitude to less than 10 -11 . In the case of HCOOH, the route via eq. [9] is assisted by 
another gas-phase route to produce HCOOH via the recombination of CH 3 Oj. This ion is 
produced following the evaporation of water via the radiative association reaction 

HCO+ + H 2 CH 3 0^ + his (11) 

although reactions of this sort have little experimental information concerning them and the 
route must be regarded as tentative at best. If this route does occur, as in our models, the 
abundance of HCOOH is not significantly affected. Thus, the abundance of HCOOH when 
the temperature reaches 100 K in both MC and RE models with E^/ Ed=0.77 is lower by ~2 
orders of magnitude in comparison with models with E^/Ej^— 0.5. Formation of the surface 
counterparts of formic acid and dimethyl ether (see Figure [T3l) occurs mainly by the accretion 
of species formed in the gas phase, and marginally by surface formation routes, which are 
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inefficient due to low mobility of radicals in models with E-^j '.Ed =0.77 '. The surface formation 
route of HCOOCH3 (eq. [5]) is inefficient in models MC-1 and MC-4, too. The low abundance 
of H 2 CO in the gas after evaporation inhibits the gas-phase formation of methyl formate via 
reaction flS]). The inefficiency of both gas phase routes and surface routes of formation leads 
to the complete absence of methyl formate in models MC-1 and MC-4 and allows only a 
small amount of the molecule in the RE model through the surface production-sublimation 
route. Below we discuss why surface formation of methyl formate via reaction (pjj) is more 
efficient in two-phase models than in multiphase models. 

As can be seen in Figures [T2l and [T3| models MC-1 and MC-4 with Eb/ED=0.3 and 
Eb/ Ed=0.5 exhibit qualitatively similar results to each other. There are two important 
differences from the results obtained with the two-phase RE model with E^/Ej^—0.5, which 
can be regarded as a reference. First, the behavior of formaldehyde and CO in the gas phase 
is more complex in multiphase models. While in the reference RE model, H2CO has a single 
desorption peak near the temperature of sublimation followed by a gradual destruction via 
gas-phase chemical reactions, in multiphase models there are two peaks. The first peak 
arises near 7 x 10 yr and corresponds to the desorption of formaldehyde from the top four 
monolayers of the grain mantle. This peak occurs with a fractional abundance of ~10 -6 . 
During the next 2 x 10 4 yr, the H 2 CO abundance drops to 10~ 9 , before the second peak 
occurs. This peak arises at ~ 1.2 x 10 5 yr and corresponds to the desorption of refractory 
water ice, which also sweeps out the remaining volatile materials. After the second peak, 
formaldehyde is gradually consumed by ion-molecule reactions with protonating ions such 
as HCO + and the neutral- neutral reaction with NH 2 , which attains a fractional abundance 
of ~ 10~ 6 after the sublimation of its precursor, ammonia, from the ice mantle. The CO 
behavior is similar, with the same two peaks in multiphase models, one of which is missing 
in two-phase models. 

The second important feature of the results obtained with multiphase models is the 
reduced abundance of gas-phase and grain surface methyl formate in comparison to the 
reference RE model. The peak abundances of both gaseous and ice methyl formate in our 
multiphase models lie up to a factor of hundred below their values in the reference RE model. 
In model MC-4 with E^/ Ed=0.5, the peak abundance of methyl formate is four times higher 
than in model MC-1. Thus, the abundance of methyl formate correlates strongly with the 
rate of surface photodissociation. The highest rate is in RE models, where all molecules in all 
~100 monolayers of mantle are available for surface photodissociation. In MC-4 models, only 
molecules in four chemically active monolayers can be photodissociated. Finally, in MC-1, 
the rate of surface photodissociation is effectively reduced by a factor of 4 in comparison 
with the MC-4 model. The rate of surface photodissociation is important for the abundance 
of HCOOCH3 because the radical precursors (OH, HCO) of surface methyl formate are 
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produced by photodissociation. 

The efficiency of the grain-surface formation route of dimethyl ether (eq. [H]) also depends 
on the number of mantle monolayers available for photodissociation, because the methyl 
radical CH3 is mainly produced by surface photodissociation of methanol. Nevertheless, the 
final gas-phase abundance of dimethyl ether is similar in two-phase and multiphase models, 
as the dominant formation route of CH3OCH3 is a gas-phase process (reaction ITU]) . 



4. Discussion 

In this study, we combined an advanced chemical model with a rather basic represen- 
tation of the evolution of a core from its prestellar to protostellar phase. It is interesting to 
compare our results with previous theoretical studies and the results of observations. 



4.1. Cold collapse phase 



Our model of the c old co llapse phase is quite similar to the model described in Section 
3.3 of iGarrod fc Paulyl (120111 ). In that study, the authors employed a modified version of a 



three- phase gas-surface-inert bulk chemical model initially proposed by lHasegawa fc Herbst 



( 19931). The var i ations of physical conditions during the cold collapse are quite similar in 
Garrod fc Paulyl (1201 ll ) and our model. We use the same density and extinction profile, but 
a somewhat differ ent profile of temperatu re. Some differences in surface chemistry are more 
significant. While IGarrod fc Paulyl (1201 if ) utiliz ed a so-called reaction -diffusion competition 
for surface reactions with activation barrie rs (IHerbst &: Millar! 120081 ) . we stick to a more 
traditional approach (lHasegawa et al.lll992l ). in which such competition is not considered. 
This difference leads to the fact that surface reactions with barriers such as CO? formation 
via the route grCO+grOH — > grC02+grH proceed more efficiently in lGarrod &: Paulyl (120 111 ) 
than in our model. Another important difference is the mechanism of oxy gen hydrogenation 
on top of a buried CO molecule implemented in IGarrod fc Paulyl ( 120111 ). We do not have 
this mechanism in our model, so the CO2 production efficiency at low temperatures in our 
model is further reduced. Nevertheless, the ice composition and structure in our stochastic 
multiphase models wi th E^/ En=0.3 and Eh/ En=0.5 are quite simi lar to those obta i ned in 
the collapse model of IGarrod fc Paulyl ( 1201 ll ). Interestingly, while IGarrod &: Paulyl (120 111 ) 
use E^/Ey)=0.3 in their collapse model, our models with the closest ice composition, MC-1, 
has E^/ En=0.5. This is because the tempera ture during collapse in our model is higher by 
a few degrees than in IGarrod fc Paulyl ( 1201 ll ). This leads to higher mobilities of species on 
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a grain surface in our model. 



The main conclusion from the comparison of our model with that of iGarrod fc Pauly 



(120 111 ) for the cold collapse phase is that the observed ice structure and composition can 
be qualitatively reproduced with an approach in which cold collapse occurs with a drop in 
temperature combined with a stochastic mult ilayer model or an approximately stochastic 
three-phase treatment (IGarrod fc Paulyl l201ll ). each of which protects already formed ice 
monolayers from further chemical processing. The adopted details of grain surface chemistry 
such as the new CO2 formation mechanism via the surface hydrogenation of oxygen on top of 
a buried CO molecule and reaction-diffusion competition are less important. The traditional 
treatment of the surface reaction grCO + grOH — > grCC-2 + H alone is sufficient to produce 
observed amounts of carbon dioxide on dust grains. 

Another control parameter of the grain surface chemistry — the diffusion/desorption 
energy ratio — is of crucial importance. As was shown above, the cold ice composition varies 
strongly in models with different adopted E^/Ed values. Unfortunately, this parameter is 
known poorly, because it is difficult to measure the diffusion barrier in laboratory experi- 
ments. It is tempting to conclude that comparison of our modeling results with observations 
of ices towards low- mass protostellar cores support a diffusion/desorption e nergy ratio of 0.5. 
Howev er, one should be careful here. First, our model as well as that of IGarrod fc Pauly 
(120 111 ) uses for comparison a me dian ice com p ositio n towards low-mass protostars from the 
infrared surveys summarized in lOberg et al.l ( 120 111 ) . While the average values are similar 
to those obtained with our model MC-1 (E h /ED=0.5), the ice composition towards indi- 
vidual sources exhibits significant variations. These variations probably reflect the different 
evolu tionary stages of individual sources and different level of energetic processing of ices 



(e.g., iGibb et al.ll2000l ). For example, there is evidence that CO ? might be more abundant 
than CO in evolved protostars such as W33a (IGibb et al.ll2000l ). Also, evolved high mass 
protostars such as W33a and GL 7009S show a signifi cantly higher than average abundance 
of methanol. The recent study of lAikawa et al.l (120121 ) also puts upper limits on the CH3OH 
abundance in the low-mass protostars L1527 and IRAS04302 to be 26% and 42% respectively 
with respect to water; these values are significantly higher than the median abundance. Such 
an ice composition is closer to our models MC-1 and MC-4 with E h / E D =0.3. Second, the ice 
composition depends on the ice thermal history during the cold collapse phase. In this study, 
a very simple temperature profile is utilized, which may differ from the thermal history of 
real ices. We did not explore this because of the huge demand of Monte Carlo models for 
CPU time. These facts indicate that comparison of our modeling results with observations is 
not sufficient to put strong constraints on the E^/En ratio. We can only claim that models 
with E^/Ey)=0.77 are eliminated by comparison with observations. 
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Another interesting result crucially dependent on t he E^/En rat i o in o ur study is the 
amount of reactive radicals trapped inside the bulk ice. iTaquet et al.l (120121 ) found in their 
model that a multilayer approach with ice porosity taken into account leads to relatively high 
abundances of reactive radicals in the bulk ice. In their reference model (constant density of 
10 5 cm -3 , constant temperature of 15 K, E^/ En=0.65), the fractional abundances of grOH, 
grHCO and grCH 3 trapped in the ice reach ~5xl0 -6 , ~10 -8 and ~10 -9 , respectively. 
We observe similarly high abundances of these radicals only in our multiphase models with 
E h /E D =0.77. Thus, in addition to the multilayer approach, which protects radicals by 
burying them in the bulk ice, a low mobility of surface species is essential, too; reactive 
radicals should be bur i ed by accreting species before they react with other species on the 
surface. ITaquet et al.l ( 120121 ) propose that abundant radicals are synthesized in the cold 



ice without photolysis, and then trapped, so that they may serve as precursors of complex 
organic species during the warm-up phase. However, as one can see at the right column of 
Fig. [12] and Fig. [131 i n our models, relatively high abundances of reactive radicals in the 
ice do not lead to the efficient formation of complex organic species in models with high 
E^/E-£,=0.77. During the warm-up, radicals are released at the temperature of sublimation 
of water ice, which is close to 100 K, and quickly evaporate where they are consumed in 
ion-molecular and photodissociation reactions. 



As shown in IVasyunin et al.l ( 120091 ) and iGarrod et al.l (120091 ) . models of surface chem- 
istry with E h /E- D =0. 3 exhibit significant stochastic effects, while models with E h / E D =0.77 
do not. Although the importance of stochastic effects in models with E^/Ed=0.5 was not 
studied in detail here, one may assume that in this intermediate case they may still be 
appreciable. If so, it means that in models that reproduce the observed ice composition, 
stochastic effects play an important role, and further modeling of ices should rely on models 
that treat them properly. Since the Monte Carlo approach utilized here is very c omputa- 
tional ly demanding, although rigorous, and the modified rate equations proposed by IGarrod 
( 120081 ) are probably not entirely correct for those aspects of surface chemistry that combine 
two-body diffusive reactions and reverse photochemical processes (see Section 2.3), the need 
for further developments in mathematical methods of astrochemical modeling is clear. 



4.2. Warm-up phase 

The MONACO code allows us for the first time to examine the impact of two features 
of ice kinetics on the chemistry in the warm-up phase leading to hot cores: the distinction 
between reactive surface vs. inert bulk and the entrapment of volatiles within water ice. We 
are particularly interested in how these feature affect the complex organic species that are 
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formed in this warm-up stage. The abundances of these species in our models are shown 
in Table [3] along with observational values. The model abundances refer to the time of 
maximum methyl formate abundance in the gas, which corresponds to 100 K. The calculated 
abundances of HCOOH, CH 3 OCH 3 and CH 3 OH are quite similar in the Monte Carlo and 
R E models with Eh/En=0 .5. The degree of agreement with observations is similar to that 



in 



Garrod fc Herbstl (120061 ) . An order of magnitude agreement is achieved for CH3OCH3 in 



all models with E^/Ed=0.5 and MC models with E^/Ed=0.3 for all sources. Agreement 
for HCOOH is reasonable for the hot corinos IRAS 16293-2422 and IRAS 4A. For the Orion 
Co mpact Ridge and Orion Hot Core, HCOOH is somewhat overproduced, as also occurs 



111 



Garrod fc Herbstl (120061 ). The abundance of gas-phase methanol after evaporation is 



higher in all models with E^/ Ed=0.5 and MC models with E^,/ Ed=0.3 by 2-3 orders of 
magnitude than observed values. Overall, the results obtained with multilayer models with 
Eb / Ej)=0 .5 and 0.3 are quite similar to the re sults obtained with our reference model RE 
with E\,/Ev=0.5, and thus similar to results of iGarrod fc Herbstl (120061 ) . 



Another picture is observed for methyl formate. While IGarrod fc Herbstl (120061 ) suc- 
ceeded in reproducing the observed abundance of HCOOCH3 with their two-phase RE model, 
our multiphase Monte Carlo models underproduce methyl formate significantly. The extent 
of the underproduction correlates with the number of mantle monolayers available for sur- 
face photodissociation. As can be seen in Figure [121 the abundance of HCOOCH3 in models 
MC-4, where four upper monolayers of the mantle are available for surface chemistry includ- 
ing photodissociation, is four times higher than in models MC-1. In two-phase RE models, 
the entire mantle population of water is available for photodissociation, which may indicate 
that photoprocessing of ices is not limited to the outermost monolayers and ta kes place at 
least in a significant fraction of the bulk ice. lAndersson fc van Dishoeckl (120081 ) found from 
theoretical calculations that photons can penetrate about a hundred monolayers and pho- 
todissociate species deeply embedded in the ice. Since the products of photodissociation 
within the bulk must somehow participate in chemical reactions, there is a need for diffusion 
within the bulk at least for some species. Bulk diffusion is not included in our current Monte 
Carlo model and is a subject for f uture research. Seve r al currently existin g models include 



a treatment of bulk diffusion (e.g., ICuppen et al.l l2009t iFayolle et al.ll201ll ). but to the best 



of our knowledge none of them has treated fully time-dependent gas-grain chemistry. 



4.3. Jumps of H 2 CO abundance in the warm-up phase 



As noted above and as can be seen in Figure [121 the abundance profile of gas-phase 
formaldehyde during warm-up is more complex in multiphase Monte Carlo models than 
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in two-phase RE models due to entrapment effects. It is intere sting to see if there is any 
observational evidence for such a complex behavior of H 2 CO. iCeccarelli et al.l ( 120011 ) re- 
ported two jumps in the gas-phase H 2 CO abundance observed around the low-mass proto- 
star I RAS16293-2422. This object is known as a "hot corino" rather than a massive "hot 
core" (IBottinelli et al.ll2004j ). However, its temperature a nd density structure are similar 
to those adopted in our model for the warm-up stage (see ICeccarelli et al.l 120001 ) . In IRAS 
16293-2422, the density rises from ~10 6 cm" 3 at 1000 AU at the edge to ~10 8 cm" 3 near 
the center of the object, while t he temperature i ncrea ses from ~40 K to 150 K over the 
same spatial range. According to ICeccarelli et al.l (120011 ). the first jump in abundance, from 
4 x 10 -10 to 4 x 10 -9 , occurs at the point where the dust temperature exceeds ~50 K, the 
sublimation temperature of formaldehyde, while the second jump, to 10 -7 , occurs at the 
point where t he dust tempera t ure re aches 100 K, the sublimation temperature of water ice. 
According to ICeccarelli et al.l (120001 ). the gas density at the first point is ~ 5 x 10 6 cm -3 , 
while at the second point it is ~ 5 x 10 7 cm -3 . For convenience, we summarize observed and 
modeled abundances of formaldehyde in Table HI From our models, we picked abundances 
of H 2 CO at three times/te mperatures, which from our point of view correspond to the three 
shells of IRAS16293-2422 (ICeccarelli et alJl200oh . Abundances in the column H 2 CO hot are 
taken from the warm-up phase of our model at time 2 x 10 5 yr and correspond to the outer 
edge of IRAS16293-2422, at a temperature of 200 K. Abundances in the column H 2 CO war m 
are taken from the warm-up phase at 1 x 10 5 yr where the temperature reaches 50 K. This 
point corresponds to the outer edge of the "warm envelope" of the source. Finally, abun- 
dances in the column H 2 CO co id are taken from the cold collapse phase of the model at the 
moment when the gas density reach es 5 x 10 6 cm" 3 . The se abundances correspond to the 
"cold envelope" of IRAS16293-2422 JCeccarelli et alJboOlh . 

One can see that the multiphase models MC-1 and MC-4 with E^/ Ed=0.3 and E^j Ed=0.5 
qualitatively reproduc e the observed H 2 CO abundance distribution for IRAS 16293-2422 
(ICeccarelli et al.l 120011 ) . Specifically, the mode led abundances of formaldehyde at cold and 
warm regions agree with ICeccarelli et al.l (120011 ) within an order of magnitude. On the other 
hand, the modeled H 2 CO abundance at the outer edge of the hot corino at 200 K is too 
high by 1-2 orders of magnitude. As compared with the multiphase mo dels, the two-phase 
RE models do not reproduce the observations of ICeccarelli et al.l (120011 ) as well, especially 
the model with E^/Ed = 0.5. The abundance of formaldehyde in the warm envelope is de- 
termined by the evaporation of H 2 CO from the ice surface while the majority of the species 
remains locked within water ice. Two-phase models cannot account for this effect. The 
exception is our RE model with E h / E D =0.3. However, this model cannot be considered as 
successful, because it produces an anomalous ice composition, as seen in Figure [TDJ Note 
that when analyzing the warm H 2 CO abundance, we ignore for two reasons the sharp spike 
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that can be seen at 8xl0 4 yr in Figure IT2l in the MC-1 and MC-4 models. First, t his spike 
occurs at a somewhat lower temperature than claimed in ICeccarelli et al.l (j200lf ) for the 
warm envelope. Second, the lifetime of this spike is very small and it should not be seen in 
observations. 



In the later study of iMaret et al.l (120041 ) . the H2CO emission was studied towards eight 
other sources. In four of them (NGC1333-IRAS 4B, NGC1333-IRAS 2, L1448-MM and 
L1527), the fitting of observational data is best when the existence of a jump in the formalde- 
hyde abundance at a temperatur e near 100 K is a s sumed . However, a second jump near 50 K, 
as found in IRAS 16293-2422 by lCeccarelli et al.l (120011 ). was not confirmed for these sources 
due to limitations in observational data. Abundances of H2CO derived from the fitting of 
observational data in IMaret et al.l ( 120041 ) for these four sources are listed in the last four 
rows of Table HI Again, the observed abundances of formaldehyde are more consistent with 
the results of multiphase than two-phase modeling. Specifically, only in multiphase models 
is the majority of H2CO released to the gas from ice mantles at temperatures near 100 K, 
as seen in Figure f~ 



5. Summary 

In this study, we constructed a multiphase gas-grain Monte Carlo model of chemical 
evolution in the interstellar medium. The model treats icy grain mantles in a layer-by-layer 
manner. The ice treatment includes distinctions between the reactive outer monolayers of 
an ice mantle and the inert bulk, as well as between zeroth and first order desorption. It 
also includes entrapment of volatile species in a refractory ice mantle. One advantage of our 
multiphase model over previously considered three-phase models is that it is capable not 
just of accreting, but also of desorbing the ice mantle in a layer-by-layer manner, which is 
important for warm-up chemistry. The behavior of ice during desorption in our model is 
justified by comparison with laboratory TPD data on ice desorption. 

With the MONACO model, the chemical evolution during the cold and warm-up phases 
leading to hot core/corino formation has been investigated. For the first time, the impact of 
a detailed multilayer approach to grain-mantle formation on the warm-up chemistry has been 
explored. Overall, our model for the first time simultaneously explains the observed compo- 
sition of ice mantles formed during the cold stage of prestellar evolution, and abundances of 
gas-phase species evaporated from the ice during the warm-up protostellar stage. 

We find that the cooling down of dust from 20 K to 10 K during the cold collapse phase, 
which is caused by the shielding of interstellar radiation during contraction from the diffuse 
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to dark cloud phase, helps us to explain the observed ice composition and structure semi- 
quantitatively. Multiphase Monte Carlo models reproduce the observed ice composition 
better than two-phase models. Also, contrary to two-phase models, multiphase models 
explain the observed ice structure, including a polar mixture of CO2 and water at monolayers 
close to the bottom of the mantle, and an apolar phase of ice rich in CO close to the top of 
the mantle. 

The impact of the multilayer ice-modeling approach on the abundances of organic 
species formed during the warm-up phase is somewhat mixed. Species such as HCOOH and 
CH3OCH3 are affected only slightly by the new approach, because their abundances to a large 
extent are controlled by gas-phase reactions. On the contrary, the abundance of HCOOCH3 
is different in two-phase and multiphase models. Here, the treatment of photochemistry is 
crucial for the modeled abundance of HCOOCH3. Unlike two-phase models, in which pho- 
todissociation occurs throughout the mantle, our multiphase Monte Carlo approach allows 
photodissociation in only the top four monolayers of the mantle. This restriction leads to an 
underproduction of gaseous methyl formate compared with observation. The failure of the 
MC models might be an indication that photochemistry and bulk diffusion could be active 
deep in the bulk ice and that some diffusion of the products to the ice surface exists in real 
ices. 

Contrary to two-phase models, o ur multiphase modeling explains jumps at 50 K and 



100 K in the H 2 CO abundance found by lCeccarelli et al.l (120011 ) in the hot corino IRAS 16293- 



2422, and jumps of H 2 CO abundance at 100 K observed in other sources by iMaret et al. 



(120041 ) . The 100 K jumps are explained as due to release of formaldehyde from the re- 
fractory water ice mantle at the sublimation temperature of H2O. This effect cannot be 
explained with two-phase models because it requires a consideration of the partial entrap- 
ment of formaldehyde in refractory water ice at 50 K. An analogous two-step behavior of 
gas-phase species in our multiphase model is predicted to occur for other volatile species 
such as CO, CH4, and CO2, and may be a subject for further observational studies of star 
formation regions. 

Various parameters controlling the grain surface chemistry have a different impact on the 
resulting ice composition. The diffusion/desorption energy ratio E^/Ed significantly affects 
the modeled ice composition. The best agreement between the modeled ice composition 
and observations is achieved in models with E^/Ej^—O.b. Models with E^/ Ej)=0.3 are also 
reasonable for the warm-up phase, while the ice composition produced with these models is 
somewhat different from observations. We included neither the competition between diffusion 
and desorption competition, nor the mechanism of ox ygen hydrogenation on top of a buried 



CO molecule, as described in iGarrod & Paulyl (1201 ll ). Nevertheless, our results are similar 
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to the results of that study for the cold ice phase. One reason for this agreement, as regards 
CO2 ice, is that the reaction grCO + grOH — > grC02 + H at the slightly higher temperatures 
we use can form CO2 on cold surfaces in sufficient amounts. 

We have used two different sets of stochastic models - MC-4 and MC-1. Models MC-1 are 
obtained by applying corrections <^ and (j3J) to model MC-4 (see Section 2.3.1). Correction 
02]) is related to the possibility of reactions between species residing in different monolayers of 
ice, while correction (EJ) is related to the deepness of penetration of photons into the ice man- 
tle. The use of models MC-1 leads to better agreement with the observed ice composition at 
the cold collapse phase than achieved with models MC-4. However, at higher temperatures, 
as occur during the warm-up stage, the two sets of models show little difference except that 
MC-4 exhibits better agreement with the observed abundance of methyl formate. 

Overall, this work shows that realistic multiphase macroscopic Monte Carlo models of 
gas-grain interstellar chemistry are not only feasible with modern computers, but also better 
explain observational data than models based on two phases and rate equations. 
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A. The MONACO algorithm 

In this section, a detailed description of the MONACO algorithm is given. Since the 



treatm ent of gas-phase chemistry in the algorithm exactly follows that described in iGillespie 



(119761 ). below we focus on the steps pertaining to surface chemistry. Nevertheless, a step- 
by-step description of the algorithm given below includes all steps, both for gas-phase and 
grain surface chemistry. 

In Figure [TH the representation of a multilayer ice mantle in the MONACO algorithm 
is shown along with assorted processes. Since a dust particle has a finite size and finite 
number of adsorption sites N s , only a finite number iV s of molecules can be placed on it in 
one monolayer. If the total number of molecules M s on a grain does not exceed N s , the 
first, or outermost, monolayer is forming. The first monolayer is considered to be chemically 
active; i.e., accretion, desorption, and chemical reactions on a grain can proceed only there. 
If M s exceeds N s , a molecule is randomly picked from the first layer and moved to an ordered 
queue with index i starting from 1. The next removed molecule will have an index i + 1 
and so on. If, due to desorption or a surface reaction with fewer products than reactants, 
the number of molecules on the outermost layer is reduced, the molecule with the highest 
index z max is moved back onto the outermost layer and i is lessened by 1. Therefore, when 
M s > N s , the total surface population of molecules in the outermost monolayer is always 
kept at N s while the number of molecules in the queue changes with time according to the 
chemical evolution of the system. In contrast to the chemically active monolayer, molecules 
in the queue are chemically inert; i.e., do not participate in chemical reactions, accretion or 
desorption processes. 

We consider the spans of z max ... z max — N s , i max — N s ... i max — 2N S etc. as the Mth 
bulk layer, (M-l)th bulk layer and so on. The layer with number 1 consists of molecules 
that came to the bulk at the beginning of system time. The last bulk layer, with the 
highest number, contains the most recent molecules buried in the bulk. Therefore, one can 
consider these model layers that together form the queue as a representation of layers in 
a physical mantle below the outermost layer. A step-by-step description of the algorithm 
is given below. The following terms are used, dj is a rate of i-th chemical reaction. S = 
(Si, 5*2, Sn, Si, Sf , S*n) is a vector of abundances of all species; Sj is an abundance of j- 
th species in gas phase, S? in an abundance of j-th species on a grain, u% is the stoichiometry 
vector for reaction k, which shows how the abundances of reactants and products of reaction 
k change due to the single event of the reaction, t current is the current physical time of a 
simulated system, t en d is the time until the chemical evolution of the system is finished, and 
N™ x is the maximum number of molecules in the chemically active monolayer. N to t is the 
actual number of molecules in chemically active monolayers. By setting N™ x to be larger 
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than N s , we assume that species in more than one physical monolayer of ice are chemically 
active. In this study, we use N™ x = 4 x N s based on the results of Section [2^2] Finally, Q and 
<imax are defined as arrays representing the ordered queue and the number of molecules in the 
queue, correspondingly. All processes in the model are treated as unimolecular or bimolecular 
reactions. In particular, accretion and desorption are represented as unimolecular chemical 
reactions: X — > grX, grX — > X. 

The algorithm steps follows: 

• Set initial abundances and the queue: S = So, Q = Qo, <?max = <?max- 

• While tcurront < ^end : 

1. Calculate rates of all reactions cij 

2. Calculate the sum of rates of all reactions: = Yll=i a i 

3. Choose uniformly distributed random numbers r x and T2 

4. Calculate time step: At = ■ ln(^-), update t curre nt = Current + ^ 

5. Choose process k to happen: k = rain : ^4=1 a « < a o • r 2 and update abundances: 
S = S + Vk 

6. If k alters surface species, update total number of molecules in the chemically 
active monolayers N tot 

7. If N tot > N™ x : 

(a) Choose uniformly distributed random number r 3 

(b) Choose molecule m to be buried in the bulk: m = rain : J2 J jZT &i < -^tot " r 3 
and update abundances: S* = - 1, g max = g max + 1, Q(q ma , x ) = m 

8. If iV to t < N™ x : return molecule r from bulk to surface: r = Q(g max ), = S* + l, 

Qmax Qmax 1 

• If tend is reached, stop the simulation and save results, else go to step [TJ 



Steps dH5] represent the classical Gillespie algorithm exactly (1 G illespid 1 1 9 761 ) . while steps 
[6HH1 are the addendum developed in this study to account for the ice composition in the 
framework of the macroscopic Monte Carlo approach. 
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Table 1: Best fit ice composition versus observed ice composition towards low-mass 
protostars a . 



Model 


Best fit time (yr) 


Min F fit 


H 2 


CO 


co 2 


CH 3 OH 


NH 3 


CH 4 


RE, E h /E D =0.3 


8.0xl0 5 


1.10 


100 


4 


170 


2 


14 


8 


RE, E h /E D =0.5 


4.0xl0 5 


0.02 


100 


41 


41 





14 


9 


RE, E h /E D =0.77 


4.7xl0 5 


0.11 


100 


49 








14 


9 


MC-1, E h /E D =0.3 


9.7xl0 5 


0.20 


100 


11 


22 


34 


25 


14 


MC-1, E h /E D =0.5 


8.3xl0 5 


0.03 


100 


42 


24 


5 


21 


9 


MC-1, E h /E D =0.77 


6.3xl0 5 


0.19 


100 


66 


2 





21 


8 


MC-4, E h /E D =0.3 


9.8xl0 5 


0.21 


100 


10 


33 


35 


24 


13 


MC-4, E h /E D =0.5 


9.8xl0 5 


0.07 


100 


30 


19 


22 


20 


10 


MC-4, E h /E D =0.77 


4.5xl0 5 


0.15 


100 


54 


4 





16 


8 


Observational value 






100 


29 


29 


3 


5 


5 



a Observational data are taken from Table 2 of lOberg et al.1 ( 2011 ) 



Table 2: Modeled peak abundances of selected gas-phase species in the cold collapse phase 
versus observational values for TMC-l(CP) a . 



Species 


Model 


Observations 


H 2 


2(-7) 


<7(-8) 


HCO+ 


4(-9) 


8(-9) 


NH 3 


3(-8) 


2(-8) 


N 2 H+ 


l(-9) 


4(-10) 


H 2 CO 


2(-8) 


5(-8) 


H 2 S 


3-9(-ll) 


5(-10) 



a Observational abundances are taken from the compilation in lSmith et al. 



(|2004h . 
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Table 3. Modeled and observed abundances of gaseous complex species in hot 

cores/corinos with respect to riH a - 



Model 


HCOOCH3-A 


HCOOCH3-E 


HCOOH 


CH3OCH3 


CH3OH 


RE, E h /E D =0.3 


2.7(-8) 




2.1(-9) 


7.3(-10) 


5.3( 


-7) 


RE, E h /E D =0.5 


2.0(-8) 




5.2(-8) 


4.3(-8) 


2.0( 


-5) 


RE, E h /E D =0.77 


7.9(-ll) 




5.8(-10) 


3.6(-12) 


4.7( 


-8) 


MC-1, E h /E D =0.3 


3.9(-10) 




4.2(-8) 


6.3(-8) 


4.1( 


-5) 


MC-1, E h /E D =0.5 


8.0(-ll) 




1.2(-8) 


1.2 (-8) 


4.8( 


-6) 


MC-4, E h /E D =0.3 


1.2(-9) 




8.9(-8) 


6.4(-8) 


3.7( 


-5) 


MC-4, £ b /£ D =0.5 


3.4(-10) 




1.4(-8) 


4.3(-8) 


2.5( 


-5) 


IRAS 16293-2422 


8.5±0.4(-8) 


1.2±0.4(-7) 


3.1(-8) 


1.2(-7) 


5.0( 


-8) 


IRAS 4A 


1.7±0.9(-8) 


1.8±0.9(-8) 


2.3(-9) 


<1.4(-8) 


<5.0( 


-9) 


Orion Compact Ridge 


1.5(-8) 




7.0(-10) 


9.5(-9) 


2.0( 


-7) 


Orion Hot Core 


7.0(-9) 




4.0(-10) 


4.0(-9) 


7.0( 


-8) 



a Models MC-1 and MC-4 with Ei>/Et)=0.77 are excluded because they produce very 
little HCOOCH3. Modeled abundances are chosen at 100 K, which is the temperature 
at which gaseous methyl formate is at its maximum abundance. The abundances of 
all species in IRAS 1 6293-2422 and IRAS 4A except for HCOOCH3-E are taken from 
Bottinelli et al.l (2007), wh ile the abundances of HCOOCH3-E in these sources are taken 



from Bottinelli et al. 



(2004J). The a bundances for the Orion Compact Ridge and Hot Core 



are taken from ISutton et al.l (119951 ) . 
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Table 4. Modeled and observed abundances of gaseous H 2 CO in hot cores/corinos with 

respect to n H a . 



Model 


hbCOhot 


H2C0 warm 


H2C0 co id 


RE, E h /E D =0.3 


4.2(-8) 


7.8(-10) 


5.2( 


-11) 


RE, E h /E D =0.5 


6.2(-7) 


1.6(-5) 


1.7( 


-11) 


RE, E h /E D =0.77 


2.6(-12) 


1.3(-8) 


3.1( 


-11) 


MC-1, E h /E D =0.3 


1.7(-6) 


3.1(-9) 


3.2( 


-11) 


MC-1, E h /E D =0.5 


9.7(-7) 


1.2(-9) 


5.0( 


-11) 


MC-1, E h /E B =0.77 


2.7(-12) 


3.7(-10) 


3.0( 


-11) 


MC-4, E h /E D =0.3 


1.6(-6) 


1.0(-8) 


4.0( 


-11) 


MC-4, E h /E D =0.5 


2.7(-6) 


4.1(-9) 


3.4( 


-11) 


MC-4, E h /E B =0.77 


4.8(-12) 


4.6(-10) 


3.7( 


-11) 


IRAS 16293-2422 


5.0(-8) 


2.0(-9) 


2.0( 


-10) 


IRAS 4B 


1.5(-6) 


2.5(-10) 






IRAS 2 


1.0(-7) 


1.5(-10) 






L1448-MM 


3.0(-7) 


4.0(-ll) 






L1527 


3.0(-6) 


1.5(-10) 







a Abundance of H 9 C O in IRAS 16293-2422 is taken from 
Ceccarelli et al.l (120011) . Abund ances for other sources are 
taken from iMaret et all (120041 ) . 
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Fig. 1. — MONACO simulation of desorption rate vs. temperature for a TPD experiment 
with 30 ML tertiary ice of composition H 2 0:C02:CO 20:1:1. Left panel: 1 monolayer is 
chemically active. Right panel: 4 monolayers are chemically active. The CO peaks are 
shown in solid black, the CO2 peaks in dashed red, and the water peak in dotted blue. 
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Fig. 2. — CO?: H?Q ratio as a meas ure of entrapment efficiency in a MONACO simulation 
of a model by iFayolle et al.l (j201l[ ) with 100 ML of a binary ice mixture of composition 
C02:H20 1:5. Left panel: 1 monolayer is chemically active. Right panel: 4 monolayers are 
chemically active. 
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Fig. 3. — Key physical parameters of the model vs. time during the cold collapse stage. 
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Fig. 4. — The evolution of the fractional abundance of simple gas-phase species during the 
whole cold collapse phase. 
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Fig. 5. — An overview of the evolution of the ice composition in different models during the 
whole cold collapse phase. The molecular concentrations are plotted as fractional abundances 
with respect to the gas-phase nuclear hydrogen density 
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t, years 



Fig. 6. — A typical plot of the calculated ice mantle thickness vs. time during the cold 
collapse phase. The grain surface is bare at the beginning of the simulation. 
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Fig. 7. — Ice evolution for the MC-1 simulation with E^/Ed— 0.3 (top row), E^/E^— 0.5 
(middle row), and E\>/ Ed=0.77 (bottom row). In the left panels, in which fraction of a 
monolayer is plotted against individual monolayers, numbered from at the surface, only 
species that contribute at least 5% of at least one monolayer at any time are shown. In 
the middle panels, where fraction of the whole mantle is plotted vs time, only species that 
contribute at least 5% of the entire bulk at any moment are shown. On histograms in the 
right column, abundances of the most abundant ice compounds with respect to solid water 
are shown at the time of minimum F^ t . 
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Fig. 8- 



The same as Figure U\ but for model MC-4. 
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Fig. 9. — The evolution of observationally important gas-phase species during the cold 
collapse phase. 
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Fig. 10. — Best fit ice composition in RE models. Left panel: Eb/ED=0.3, middle panel: 
E b /E D =0.5, right panel: E b /E D =0.77. 
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Fig. 11. — Abundances of major gas-phase species during the warm-up phase. 
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Fig. 12. — Abundances of organic gas-phase species during the warm-up phase. 
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Fig. 13. — Abundances of organic species on the grain surface during the warm-up phase. 
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Fig. 14. — Processes occurring on an ice mantle in which the outermost layer is full in the 
MONACO model. Excess species adsorbed into this layer or produced on it by reaction lead 
to the random addition of an excess species to a virtual queue, shown beneath the outermost 
layer. The virtual queue is reduced in its number of species when desorption occurs from 
the outermost layer. 



